/*!
 ***************************************************************************
 * \file   heat3.cpp 
 * Solving the Heat difusion equation in 3D.
 ***************************************************************************
 *  \author    Luis M. de la Cruz Salas [ Fri Apr 18 08:58:33 BST 2008 ]
 ***************************************************************************
 * This code solves the heat diffusion equation in three dimensions: 
 * \f$ \frac{\partial T}{\partial t} = \nabla^2 T\f$. 
 * \par
 * This equation is solved in a unit cubic prism \f$ x, y, z \in [0,1] \f$. 
 * The boundary conditions are as shown in the next figure:
 *
 * \image html cubo3D.png "Unit cubic prism" width=10cm
 * \image latex cubo3D.eps "Unit cubic prism" width=10cm
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./heat3  \endverbatim
 * The results are stored in \c Data3 .
 ***************************************************************************
 */
#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/CDS.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"

using namespace Tuna;

// Scalar field in three dimensions
typedef TunaArray<double,3>::huge ScalarField3D;

int main() 
{
    double length_x, length_y, length_z;
    int num_nodes_x, num_nodes_y, num_nodes_z;
    int num_vols_x, num_vols_y, num_vols_z;
    double tolerance, dt;
    int frequency, tdma_iter, max_iter;

// Reading data from the input file
    std::ifstream input_cin ("input3"); 
    input_cin >> length_x;
    input_cin >> length_y;
    input_cin >> length_z;
    input_cin >> num_nodes_x; 
    input_cin >> num_nodes_y; 
    input_cin >> num_nodes_z; 
    input_cin >> dt;
    input_cin >> max_iter;
    input_cin >> tolerance;
    input_cin >> tdma_iter;
    input_cin >> frequency;
    input_cin.close();

    StructuredMesh<Uniform<double, 3> > mesh(length_x, num_nodes_x,
					     length_y, num_nodes_y,
					     length_z, num_nodes_z);
    float dx = mesh.getDelta(X);
    float dy = mesh.getDelta(Y);
    float dz = mesh.getDelta(Z);
    mesh.print();

    ScalarField3D T(mesh.getExtentVolumes());
    ScalarField3D Tn(mesh.getExtentNodes());
    
// Boundary condition : Dirichlet
    Range all = Range::all();
    T(all, T.lbound(secondDim), all) = 0.5;  // Bottom wall
    T(all, T.ubound(secondDim), all) = -0.5; // Top wall

    NumUtils::interpolateToNodes(Tn, T);

#ifdef WITH_DX
    InOut::writeToFile_DX(Tn, 0, "./Data3/temp.", dx, dy, dz);
#endif

#ifdef WITH_AM
    InOut::writeToFile_AM(Tn, 0, "./Data3/temp.", dx, dy, dz);
#endif

    SparseMatrix< Diagonal<double, 3> > A(num_nodes_x, num_nodes_y, num_nodes_z); 

    ScalarField3D               b(num_nodes_x, num_nodes_y, num_nodes_z);

    ScalarEquation<CDS<double, 3 > > heat(T, A, b, mesh.getDeltas());
    heat.setDeltaTime(dt);
    heat.setGamma(1.0);
    heat.setDirichlet(TOP_WALL);
    heat.setDirichlet(BOTTOM_WALL);
    heat.setNeumann(LEFT_WALL);
    heat.setNeumann(RIGHT_WALL);
    heat.setNeumann(FRONT_WALL);
    heat.setNeumann(BACK_WALL);
    heat.print();

    ofstream error_fp("./Data3/error"), residual_fp("./Data3/residual");
    int iteration = 0;
    double error = 1.0, residual = 1.0;
    
    while ( (error > tolerance) && (iteration < max_iter) )
    {
	heat.calcCoefficients();
	Solver::TDMA3D(heat, tolerance, tdma_iter);
	error = NumUtils::calcErrorL2(T, heat.phi);
	residual = heat.calcResidual();

	error_fp << iteration << "\t" << error << endl;
	residual_fp << iteration << "\t" << residual << endl;

	heat.update();

	cout << "\n ---- Iteration = " << ++iteration
	     << " Residual = " << residual
	     << " Error = " << error;
	if ( !(iteration % frequency) ) {
	    NumUtils::interpolateToNodes(Tn, T);

#ifdef WITH_DX
	    InOut::writeToFile_DX(Tn, iteration, "./Data3/temp.", dx, dy, dz);
#endif

#ifdef WITH_AM
	    InOut::writeToFile_AM(Tn, iteration, "./Data3/temp.", dx, dy, dz);
#endif
	}
    }

#ifdef WITH_DX
    std::cout << "\n\n NOTE: The data were written in the [./Data3] directory in DX format";
#endif
    return 0;
}




